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Abstract 



O 

in 
o 

■ We study the magnetic interaction between a neutron star and its disk by solving 

Oh 1 the time-dependent relativistic force-free equations. At the initial state, we assume 

that the dipole magnetic field of the neutron star connects the neutron star and 
^ \ its equatorial disk, which deeply enters into the magnetosphere of the neutron star. 

Magnetic fields are assumed to be frozen to the star and the disk. The rotation of 



the neutron star and the disk is imposed as boundary conditions. We apply Harten- 
Lax-van Leer (HLL) method to simulate the evolution of the star-disk system. We 
carry out simulations for (1) a disk inside the corotation radius, in which the disk 
rotates faster than the star, and (2) a disk outside the corotation radius, in which the 
neutron star rotates faster than the disk. Numerical results indicate that for both 
models, the magnetic field lines connecting the disk and the star inflate as they are 
twisted by the differential rotation between the disk and the star. When the twist 
angle exceeds 7r radian, the magnetic field lines expand with speed close to the light 
speed. This mechanism can be the origin of relativistic outflows observed in binaries 
containing a neutron star. 

Key words: accretion, accretion disks — method: numerical — stars: magnetic 
fields — stars: winds, outflows 
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1. Introduction 



Magnetic fields play essential roles in various activities observed in radio pulsars and in 
X-ray pulsars. In X-ray pulsars, an accretion disk surrounding a neutron star interacts with the 
magnetic field of the neutron star. In low mass X-ray binaries (LMXBs), angular momentum 
supplied from the accretion disk accelerates the rotation of the neutron star. These objects are 
sometimes observed as rejuvenated radio pulsars (e.g., Backer et al. 1982). Even in isolated 
radio pulsars, recent observations by the Chandra X-ray satellite (Weisskopf et al. 2000; Hester 
et al. 2002; Pavlov et al. 2003) revealed the existence of disk-like distribution of hot plasmas 
around the pulsar and the existence of bipolar jets. Magnetic interaction between the neutron 
star and its surrounding disk can be the origin of jet-like outflows. 

The interaction between a magnetized star (neutron star, white dwarf or a young stellar 
object) and its accretion disk has been studied extensively. When magnetic field lines connect 
the star and its disk, their torque spins up or down the central star (e.g., Ghosh, Lamb 1978, 
1979a, 1979b; Lovelace et al. 1995). 

Zylstra (1988) obtained a sequence of twisted force- free fields around a neutron star and 
showed that when critical twist is accumulated, loss of equilibrium leads to the expansion of the 
magnetosphere. Lynden-Bell and Boily (1994) used semianalytic techniques for non-relativistic 
force-free configurations. By studying the evolution of force-free magnetic loops anchored to 
the star and the disk, they obtained self-similar solutions for it. They found that twist injection 
from the rotating disk makes the magnetic loops unstable and inflate and that the loops expand 
along the direction 60° from the rotational axis of the disk. Lovelace et al. (1995) applied this 
mechanism to keplerian disks. 

In non-relativistic MHD case, Hayashi, Shibata, and Matsumoto (1996) carried out 
two-dimensional (axisymmetric) magnetohydrodynamic simulations of magnetic interaction be- 
tween a star and its disk. They showed that magnetic reconnection takes place in the current 
sheet created inside the inflating magnetic loops. Their model can explain the X-ray flares 
and outflows observed in protostars. Miller and Stone (1997) presented the results of simula- 
tion including the rotation of the central star. Goodson et al. (1997, 1999) and Goodson and 
Winglee (1999) carried out longer time scale simulations of star-disk magnetic interaction and 
demonstrated quasi-periodic formation and destruction of magnetosphere by twist injection 
from the disk. Kato, Hayashi, and Matsumoto (2004) reported the results of MHD simulation 
of the magnetic interaction between a neutron star and its disk. They showed that due to the 
gas pressure external to the expanding loop, magnetic tower (Lynden-Bell 1996) is formed and 
that semi-relativistic jet flows out along the magnetic tower. This mechanism can explain the 
origin of semi-relativistic jets observed in Sco X-l (Fomalont et al. 2001a, 2001b). In their 
simulation, however, special relativistic effects were neglected. 

In magnetospheres of neutron stars, relativistic effects are not negligible. The purpose of 
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this paper is to extend the disk-star magnetic interaction model to a relativistic regime. When 
plasma density is low or magnetic field is strong, the Alfven speed can be close to the light speed. 
Thus, we need to solve relativistic equations to study the evolution of the magnetosphere. Here 
we present numerical results for the magnetosphere dominated by the electromagnetic field, 
which can be described by time dependent relativistic force-free equations. 

Recently, Komissarov (2002) carried out simulations of relativistic force-free fields to 
study the electrodynamics of pulsar magnetospheres and black holes. Here we extend the 
model for the case including a rotating equatorial disk threading the magnetosphere. 

In section 2, we describe the force-free equations, numerical scheme, and redefinition 
of velocity fields. The results of simulations are given in section 3. Section 4 is devoted for 
summary and discussion. 

2. Numerical models 

2.1. Basic equations 

We use force-free equations (e.g. Uchida 1997 for basic theory of force-free fields) derived 
from the special relativistic magnetohydrodynamics (SRMHD) equations to simulate the pulsar 
magnetosphere. The force-free approximation is applicable when the electromagnetic energy 
density is much larger than the energy density of the plasma. The force-free equations are 

f + V.M = 0, (1) 

1 »D 

+ Vx£ = 0, (2) 



where 



c dt 

1 



Aire 

is momentum, 

47T 



E x B (3) 



E i E j + B i B j - -S ij (E-E + B-B) 
2 v J 



(4) 



is electromagnetic stress tensor and c is speed of light. Equations (1), (2), (3) and (4) are 
equivalent to Maxwell's equations and the force-free condition, 

p e E + - J x B = 0, (5) 

c 

where p e is charge density and J is current density. Equation (5) implies degeneracy of electr- 
magnetic fields, 

E-B = 0, (6) 

which replaces the perfect conductivity condition in classical MHD (see Komissarov 2002). 
Under this condition, we can solve equation (3) as 



3 



We adopt P and B as fundamental variables for time dependent simulation. We normalize 
the equations by using the typical strength of magnetic field Bq (on the surface of central star 
at rotational axis), and radius of the central star Rq as units. The unit of time is tq = Rq/c. 
We use a spherical coordinate system (r,6,<p) with 6 = parallel to the star's rotational axis. 
We assume axisymmetry (d/d<f) = 0) but all three components of P and B are retained. We 
assume that the star's dipole moment is aligned with the rotational axis. 

2.2. Definition of velocity fields 

We introduce the velocity of magnetic field lines «g as 



This determines the velocity of the magnetic field lines uniquely. By using (8), the electric field 
E is written as 

E=--v B xB. (9) 

When the electric field is represented like equation (9), the electromagnetic field must satisfy 

B 2 -E 2 >0. (10) 

We checked that the condition (10) is satisfied throughout the calculation. Since the force-free 
equations are derived from the relativistic MHD equations by retaining the electromagnetic field 
only, plasma's velocity v does not appear explicitly in the theory. However, we can assume the 
existence of the plasma that is frozen to the magnetic fields. If equation (10) is satisfied, the 
plasma can flow along the magnetic field lines. 

2.3. Numerical scheme 

We use the HLL method (Harten et al. 1983; Janhunen 2000) to solve the time dependent 
equations. It belongs to the family of an upwind scheme. In contrast to the Roe type scheme 
which requires eigenvalues and eigenvectors of the characteristic matrix, HLL scheme only needs 
eigenvalues. We used MUSCL-type method to attain second order accuracy in space. The HLL 
method is formulated as follows (Harten et al. 1983). A one dimensional hyperbolic system 
can be written in the form 

du 9F 3, (11) 



dt dx 

where the vector of conservative variables is 

U=(P\B l Y, (12) 

F = F(U) is the flux vector, and S = S(x,t) is the source term. The fluxes are defined at 
cell interfaces. We denote the interface variables at the right- and left-hand side of each cell 
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interface as E/ R and Ul. We denote the maximum speed of right- and left-going wave as & R and 

6l (in this simulation, 6 R = +c and 6l = — c) and the fluxes as .F R = F(Ur) and .Fl = F([/l). 

The HLL flux is then given by 

6 L F R + 6 R F L -6 R 6 L (t/ R -C/ L ) 

-thll— ; ; • {i-o) 

Or ~ b L 

Second order accuracy in time is achieved by first advancing U by At/ 2 and next advancing 
U by At using the flux at At/2. 

We use 450 (model 1) and 460 (model 2) non-uniform grids in the radial direction 
and 180 uniform grids in ^-direction. The grid size is Ar = 0.0125 when r < 1.125i?o and 
gradually increases as Ark+i/ Ark = 101 when r > 1.125i?o- The maximum radius is r max = 
100i?o (model 1) and r max = 110.Ro (model 2). 

2.4- Boundary conditions 

We assume that dipole magnetic fields initially connect the central star and the geomet- 
rically thin disk at the equatorial plane (3.0 < r < 95.0). We assume that the central star and 
the disk are perfect conductors. At the equatorial plane, we impose boundary conditions such 
that B r and B$ are antisymmetric and Bg is symmetric with respect to 9 = tt/2. For electric 
fields, we impose the following conditions at 9 = n/2 (3.0 < r < 95.0), 

E r = — Be, 

E e = --B r , (14) 

c 

^ = 0, 

which imply that magnetic fields rotate with the disk. Here u is the angular velocity of the 



Keplerian disk (tu = ^GM/r 3 , where G and M are gravitational constant and mass of the 
central star, respectivery). At the equatorial plane in 1.0 < r < 3.0, we imposed the condition 
that E r and are symmetric and Eg is antisymmetric with respect to 9 = tt/2. The magnetic 
field component Bg at the equatorial plane is extraporated from Bg at the grid point next to 
the equator. The electric field E r at the equator is computed by using equation (14). The 
momentum P at the equator is computed by using equation (3) and (14). 
The electric field on the stellar surface satisfies 
R ttsin9 

tj r = Dg, 

E = flossing (15) 

c 

^ = 0, 

where Q is the angular velocity of the central star. The magnetic fields B r , Bg and B^ at the 
stellar surface are extraporated from B r , Bg and B^ at the point next to the stellar surface. The 
electric field at the stellar surface is computed by using equation (15). The momentum P is 
computed by using equation (3) and (15). In order to check that the footpoints of magnetic field 
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Fig. 1. Time evolution of magnetic field lines projected onto the r — 6 plane (solid curves) in model 1. 
The color scale shows the twist B,p/B p . Arrows show the poloidal component of the velocity of magnetic 
field lines. 

lines don't slip on the surface of the star and the disk, we computed the velocity of magnetic 
field lines near the surface of the star and the disk by using equation (8). We confirmed that 
vb coincides well with the rotational velocity of the accretion disk and the stellar surface. 

We carried out simulations for two models. In model 1, we take Q = and uo = 
0.3(c/Ro)(r/R )~ 3/2 . In model 2, we take Q = 0.05c/ Rq and u = 0. Model 2 corresponds 
to the disk far outside the corotational radius. We use free boundary conditions at the outer 
boundaries (r = r max ). On the rotational axis (8 = 0), we substituted the values P r and B r at 
the grid point next to the rotational axis. 

3. Numerical results 

Results of simulations for model 1 are shown in figure 1. Solid curves in figure 1 show 
the magnetic field lines projected onto the r-9 plane. The color scale shows the strength of 
twist of magnetic field (B^/B p ) where B^ is the toroidal magnetic field and B p is the poloidal 
magnetic field. The arrows show the poloidal component of the velocity of magnetic field lines 
defined by equation (8). 

The disk at r = 3Rq rotates about 7 radian at t = 120.0. By the rotation of the disk, 
the magnetic field lines connecting the central star and the disk are twisted. Torsional Alfven 
waves propagating along magnetic field lines are reflected on the surface of the central star. 
The typical crossing time of Alfven waves travelling between the central star and the disk is 
~ 2r/c. As time goes on, oscillations driven by the torsional Alfven wave dissipates. Since the 
toroidal component of magnetic fields increases near the central star by accumulation of twist, 
magnetic pressure increases. From about t = 50.0, magnetic field lines begin to inflate. This 
rapid expansion of twisted magnetic loops takes place when the twist angle exceeds about it 
radian. The inflation is due to the loss of equilibrium in the twisted magnetic loops. Strong 
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Fig. 2. The same as figure 1 but for model 2. Solid curves show magnetic field lines. Color scale shows 
the twist B^/B v . Arrows show the poloidal component of the velocity of magnetic field lines. 
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Fig. 3. (Left) Time evolution of the maximum distance from r = to the point on the magnetic field line 
anchored to the disk at r = 4.9. The solid curve is for model 1 and the dotted curve is for model 2. The 
dashed line represents the wave propagating at light speed. (Right) Time evolution of the twist B^/Bp at 
the same point as that in the left figure. The solid curve is for model 1 and the dotted curve is for model 
2. 

current layer is created inside the expanding magnetic loops. These results are consistent with 
the stability analysis by van Ballegooijen (1994) and Lynden-Bell and Boily (1994). Figure 2 
shows the time evolution of magnetic field lines (solid curves), twist (B^/ B p ), and the poloidal 
component of the velocity of magnetic field lines (arrows) for model 2. In this model, the central 
star is rotating faster than the equatorial disk. Similarly to the model 1, the magnetic loops 
rapidly expand after t = 50 and form a long thin current sheet. The solid and dotted curves 
in figure 3 (left) show the maximum distance from r = to the magnetic field line anchored to 
the equatorial disk at r = 4.9. The expansion speed approaches the light speed (dashed line). 
Figure 3 (right) shows the degree of twist (B^/Bp) at the same point as that of figure 3 (left). 
The twist increases with time until t ~ 50. Subsequently, the magnetic field lines inflate almost 
keeping the degree of twist. Figure 4 shows the poloidal velocity of the magnetic field line at 
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Fig. 4. The poloidal velocity of the magnetic field line at the same point as those depicted in figure 3. 
The velocity of magnetic field lines is defined by equation (8). The solid curve is for model 1 and the 
dotted curve is for model 2. 
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Fig. 5. The dependence of the radial expansion of magnetic fields on the angular speed of the central 
star, fi, when ui = in the model 2. The curves show the maximum distance from r = to the point on 
the magnetic field line anchored to the disk at r = 4.9. The solid curve is for model 2. The dotted curve 
is for fl = 0.075, and the dashed curve is for = 0.025. 

the same point as those depicted in figure 3. The poloidal velocities in model 1 and 2 exceed 
90% of the light speed at t ~ 100 and gradually approach the speed of light. Consequently, the 
toroidal velocity decreases and approach to zero. 

Figure 5 shows the dependence of the radial expansion of magnetic fields on the angular 
speed of the central star, Q, when uo = 0. The solid curve is for model 2. As Q increases, the 
poloidal component of the velocity of magnetic field lines approaches the light speed within 
shorter period. 

4. Summary and Discussion 

In this paper, we numerically solved the relativistic force-free degenerate electrodynamic 
equations to study the time evolution of magnetic field lines connecting a neutron star and its 
disk. This approach is an extension of the study by Uzdensky et al. (2002) in which they solved 
force-free magnetosphere by using semianalytic methods. Starting from potential magnetic 
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fields, we showed that magnetic fields are twisted by the differential rotation between the disk 
and the central star. After critical twist accumulates, the magnetic field lines connecting the 
star and the disk inflate with speed exceeding 0.95c. In the region distant from the central star 
and the disk, the poloidal velocity approaches the speed of light. When the plasma is loaded to 
these expanding magnetic loops, we can expect relativistic outflows with Lorentz factor T > 3. 
A long thin current sheet is created inside these expanding magnetic loops. The elongation of 
the current sheet continues until the end of the calculation. 

In this paper, we have solved ideal force-free equations. When the resistivity is included, 
magnetic reconnection will occur in the current sheet. Since force-free equations can't handle 
magnetic reconnection, we need to solve the full set of resistive MHD equations to simulate 
the magnetic reconnection. Another limitation of force-free approximation is that gas pressure 
does not appear in force-free equations. The existence of plasma pressure or magnetic fields 
helps collimating the expanding magnetic field lines toward the rotational axis (e.g. Kato et al. 
2004; Tsinganos, Bogovalov 2002). Relativistic MHD simulations including gas pressure will be 
able to reproduce the collimated relativistic jets. We will report the results of such simulations 
in subsequent papers. 

The authors thank Y. Uchida, S. Hirose, M. Nakamura for discussion. We thank 
S. Komissarov for discussion during the visit of R. Matsumoto to the University of Leeds. 
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by the grant of ministry of education, science, sports, culture and technology (15037202, P.I. 
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